{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Navier Stokes Equations\n",
    "===\n",
    "Find velocity $u : \\Omega \\times [0,T] \\rightarrow R^d$ and pressure $p : \\Omega \\times [0,T] \\rightarrow R$ such that\n",
    "\n",
    "$$ \n",
    "\\begin{array}{ccccl}\n",
    "\\frac{\\partial u}{\\partial t} - \\nu \\Delta u + u \\nabla u & + & \\nabla p & = & f \\\\\n",
    "\\operatorname{div} u & & & = & 0\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import SplineGeometry\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schäfer-Turek benchmark geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (2, 0.41), bcs = (\"wall\", \"outlet\", \"wall\", \"inlet\"))\n",
    "geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc=\"cyl\", maxh=0.02)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.07))\n",
    "\n",
    "mesh.Curve(3)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Higher order Taylor-Hood element pairing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = VectorH1(mesh,order=3, dirichlet=\"wall|cyl|inlet\")\n",
    "Q = H1(mesh,order=2)\n",
    "X = FESpace([V,Q])\n",
    "\n",
    "u,p = X.TrialFunction()\n",
    "v,q = X.TestFunction()\n",
    "\n",
    "nu = 0.001  # viscosity\n",
    "stokes = nu*InnerProduct(grad(u), grad(v))+ \\\n",
    "    div(u)*q+div(v)*p - 1e-10*p*q\n",
    "\n",
    "a = BilinearForm(X)\n",
    "a += stokes*dx\n",
    "a.Assemble()\n",
    "\n",
    "# nothing here ...\n",
    "f = LinearForm(X)   \n",
    "f.Assemble()\n",
    "\n",
    "# gridfunction for the solution\n",
    "gfu = GridFunction(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "parabolic inflow at inlet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "# Draw (Norm(gfu.components[0]), mesh, \"velocity\", sd=3)\n",
    "Draw (gfu.components[0], mesh, \"vel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "solve Stokes problem for initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv_stokes = a.mat.Inverse(X.FreeDofs())\n",
    "\n",
    "res = f.vec.CreateVector()\n",
    "res.data = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += inv_stokes * res\n",
    "\n",
    "scene = Draw (gfu.components[0], mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "implicit/explicit time-stepping:\n",
    "\n",
    "$$\n",
    "\\frac{u_{n+1}-u_n}{\\tau} - \\nu \\Delta u_{n+1} + \\nabla p_{n+1} = f - u_n \\nabla u_n\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\operatorname{div} u_{n+1} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = 0.001 # timestep parameter\n",
    "\n",
    "mstar = BilinearForm(X)\n",
    "mstar += (u*v+tau*stokes)*dx\n",
    "mstar.Assemble()\n",
    "inv = mstar.mat.Inverse(X.FreeDofs(), inverse=\"sparsecholesky\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the non-linear convective term $\\int u \\nabla u v$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conv = BilinearForm(X, nonassemble = True)\n",
    "conv += (Grad(u) * u) * v * dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "implicit Euler/explicit Euler splitting method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 0\n",
    "tend = 1\n",
    "i = 0\n",
    "gfut = GridFunction(V, multidim=0)\n",
    "vel = gfu.components[0]\n",
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        # print (\"t=\", t, end=\"\\r\")\n",
    "\n",
    "        conv.Apply (gfu.vec, res)\n",
    "        res.data += a.mat*gfu.vec\n",
    "        gfu.vec.data -= tau * inv * res    \n",
    "\n",
    "        t = t + tau\n",
    "        i = i + 1\n",
    "        if i%10 == 0:\n",
    "            scene.Redraw()\n",
    "            gfut.AddMultiDimComponent(vel.vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Draw (gfut, mesh, interpolate_multidim=True, animate=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
